function H=load3DHFilter(filename,precision)
    global verbose
    fid=fopen(filename,'r');
    if fid > 0
        n=fread(fid,1,'int');
        timestamp=fread(fid,1,'real*8');
        precond=fread(fid,1,'int');
        if verbose
            fprintf('loading 3D filter which was created on %s for n = %d and precond %d\n',datestr(timestamp),n,precond);
        end
        if n > 2
            K=n+2*nchoosek(n,2)+nchoosek(n,3);
        else
            K=n+2*nchoosek(n,2);
        end
        if strcmp(precision,'single')
            Htmp0=fread(fid,K,'real*4');
        elseif strcmp(precision,'double')
            Htmp0=fread(fid,K,'real*8');
        else
            error('precision must be single or double');
        end
        fclose(fid);
        Htmp1=zeros(n,n,n,precision);
        
        W=calculateIndexMapping(n,precision);
        k=1;
        for i = 1:K
            j=1;
            while j <= 6 && W(i,j) > 0
                Htmp1(W(i,j))=Htmp0(k);
                j=j+1;
            end
            k=k+1;
        end
        Htmp2=cat(1,Htmp1,flipdim(Htmp1(1:n-1,1:n,1:n),1));
        Htmp3=cat(2,Htmp2,flipdim(Htmp2(1:2*n-1,1:n-1,1:n),2));
        H.filter=cat(3,Htmp3,flipdim(Htmp3(1:2*n-1,1:2*n-1,1:n-1),3));
        H.precond=precond;
    else
        error('failed to open file %s for reading, filter was not loaded\n',filename);
    end
end

function W=calculateIndexMapping(n,precision)
    if n > 2
        V=zeros(nchoosek(n,3),5,precision);
    end
    U=zeros(2*nchoosek(n,2),2,precision);
    if n > 2
        W=zeros(n+2*nchoosek(n,2)+nchoosek(n,3),6,precision);
    else
        W=zeros(n+2*nchoosek(n,2),6,precision);
    end
    if n > 2
        t=1;
        for i=1:n-2
            k=1;
            for j=i:-1:1
                s=k;
                for m=1:k
                    V(t,:)=[j,n*s-j,s+j,n*j-s,s];
                    t=t+1;
                    s=s-1;
                end
                k=k+1;
            end
        end
    end
    t=1;
    for i=1:n-1
        for j=i:-1:1
            U(t,:)=[j,n*j];
            t=t+1;
        end
        for j=i:-1:1
            U(t,:)=[n*j,j];
            t=t+1;
        end
    end
    if n > 2
        V=(n-1)*V;
    end
    U=(n-1)*U;
    p=1;
    for i=1:n
        for j=0:i-1
            for k=0:j
                W(p,1)=i+j*n+k*n^2;
                p=p+1;
            end
        end
    end
    p=1;
    t=1;
    s=1;
    for i=1:n-1
        p=p+1;
        for q=1:2
            W(p,q+1)=W(p,q)+U(t,q);
        end
        p=p+1;
        t=t+1;
        for k=1:i-1
            if n > 2
                for m=1:k
                    for q=1:5
                        W(p,q+1)=W(p,q)+V(s,q);
                    end
                    p=p+1;
                    s=s+1;
                end
            end
            for q=1:2
                W(p,q+1)=W(p,q)+U(t,q);
            end
            p=p+1;
            t=t+1;
        end
        for q=1:2
            W(p,q+1)=W(p,q)+U(t,q);
        end
        p=p+1;
        t=t+1;
        for j=i-1:-1:1
            for q=1:2
                W(p,q+1)=W(p,q)+U(t,q);
            end
            p=p+1;
            t=t+1;
        end
    end
end